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The superorbital variability and triple nature of the X-ray source 
4U 1820-303 
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ABSTRACT 

We perform a comprehensive analysis of the superorbital modulation in the ultracompact X- 
ray source 4U 1820-303, consisting of a white dwarf accreting onto a neutron star. Based on 
RXTE data, we measure the fractional amplitude of the source superorbital variability (with 
a ^170-d quasi-period) in the folded and averaged light curves, and find it to be by a factor 
of ~2. As proposed before, the superorbital variability can be explained by oscillations of the 
binary eccentricity. We now present detailed calculations of the eccentricity-dependent flow 
through the inner Lagrangian point, and find a maximum of the eccentricity of ~ 0.004 is suf- 
ficient to explain the observed fractional amplitude. We then study hierarchical triple models 
yielding the required quasi-periodic eccentricity oscillations through the Kozai process. We 
find the resulting theoretical light curves to match well the observed ones. We constrain the 
ratio of the semimajor axes of the outer and inner systems, the component masses, and the 
inclination angle between the inner and outer orbits. Last but not least, we discover a remark- 
able and puzzling synchronization between the observed period of the superorbital variability 
(equal to the period of the eccentricity oscillations in our model) and the period of the general- 
relativistic periastron precession of the binary. 

Key words: accretion, accretion discs - binaries: general - globular clusters: individual: NGC 
6624 - stars: individual: 4U 1820-303 - X-rays: binaries - X-rays: stars. 



1 INTRODUCTION 

4U 1820-303 is one of the most remarkable low-mass X-ray bi- 
naries. This ultracompact binary consists of an A/2 = (0.06- 
0.08)Mq He white dwarf secondary (Rappaport et al. 1987, here- 
after R87) accreting via Roche-lobe overflow onto a neutron star 
(of the mass Mi). Its binary period is as short as Pi ~ 685 s 
(Stella, Priedhorsky & White 1987; Anderson et al. 1997). A likely 
scenario for the formation of 4U 1820-303 appears to be a direct 
collision of a neutron star and a giant (Verbunt 1987; Ivanova et 
al. 2005) in its parent globular cluster NGC 6624. This resulted in 
a binary consisting of the neutron star and the stripped giant core, 
which then cooled down to the present relatively degenerate state 
(R87). The secondary has to have a very low H abundance in or- 
der to fit its Roche lobe (R87), and the occurence of type-I X-ray 
bursts implies that it cannot be made of elements heavier than He. 
Interestingly, 4U 1820-303 was the first X-ray burster identified 
with a known X-ray source (Grindlay et al. 1976). The most likely 
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distance to the source appears to be 7.6 ± 0.4 kpc (Kuulkers et al. 
2003), as discussed in a companion paper (Zdziarski et al. 2007, 
hereafter Z07). 

The accretion occurs as a consequence of the loss of the angu- 
lar momentum of the binary via emission of gravitational radiation 
(Paczynski 1967), as proposed for this system by Stella et al. (1987) 
and Verbunt (1987), and calculated in detail by R87. The implied 
present mass-loss rate from the secondary (R87; see Section[3]be- 
low) is completely compatible with that corresponding to the aver- 
age isotropic bolometric luminosity measured from the Rossi X-ray 
Timing Explorer (RXTE) data by Z07. This provides a strong sup- 
port for this accretion model. The present phase of the high mass 
transfer started only ~ 10 6 yrs ago and may be sustained only for 
~3 x 10 B yrs (R87). 

The 685-s period was discovered in X-rays as a modulation 
with a ~2-3 percent peak-to-peak amplitude (Stella et al. 1987). 
The period is very stable, with a low Pi jP\ = (— 3.5±1.5) x 10~ 8 
yr _1 (Chou & Grindlay 2001, hereafter CG01), which makes it 
certain it is due to the binary motion. The modulation was proposed 
to be due to the obscuration of the X-ray source by a structure at 
the edge of the accretion disk. A stronger modulation in the UV was 
predicted by Arons & King (1997), and subsequently discovered by 
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Anderson et al. (1997), with Pi = 687.6 ± 2.4 s and a 16 per cent 
peak-to-peak; amplitude. 

The most unusual feature of 4U 1820-303 appears to be the 
luminosity variation by factor of > 2 at a (superorbital) period 
of P 3 ~ 170 d (Priedhorsky & Terrell 1984; Smale & Lochner 
1992; CG01; Simon 2003; Wen et al. 2006, hereafter W06). CG01 
found the modulation is stable with P3 = 171.0 ± 0.3 d and 
|P 3 /P 3 | < 2.2 xl0 _4 yr _1 based on data from 1969 to 2000. W06 
found P 3 = 172 ± 1 d based on 8.5 yrs of RXTE All Sky Moni- 
tor (ASM) data. The fact that X-ray bursts take place only at the 
flux minima (Cornelisse et al. 2003, Z07) proves that the observed 
variability is due to intrinsic luminosity/accretion rate changes and 
not due to, e.g., obscuration or changes of the projected area of the 
source due to precession. This is further supported by strong cor- 
relations between the observed flux and the source spectral state, 
varying with the flux in a way typical of atoll-type neutron-star bi- 
naries (Bloser et al. 2000; Gladstone, Done & Gierlihski 2007), and 
between the frequency of kHz QPOs observed from the source and 
the flux (Zhang et al. 1998; van der Klis 2000). The source was 
classified as an atoll by Hasinger & van der Klis (1989). 

The intrinsic, accretion-rate related, character of the long-term 
periodic flux changes in 4U 1820-303 is unique among all sources 
showing superorbital variability. In other cases, such changes of 
the observed flux appear compatible with being caused by accre- 
tion disc and/or jet precession, which either results in variable ob- 
scuration of emitted X-rays as in Her X-l (Katz 1973), or changes 
the viewing angle of the presumed anisotropic emitter, as in SS 
433 (Katz 1980) or Cyg X-l (e.g., Lachowicz et al. 2006), or both. 
However, that precession keeps the inclination angle of the disc 
with respect to the orbital plane constant, and thus it cannot change 
the accretion rate (or the luminosity). Also, the ratio between the 
superorbital and orbital periods is ~ 2.2 x 10 4 , which is much 
higher than that possible to obtain from any kind of disc preces- 
sion at the mass ratio of the system (e.g., Larwood 1998; Wijers & 
Pringle 1999). This appears to rule out also models which would 
attempt to explain the variable accretion rate by a varying disc in- 
clination angle (with respect to the orbital plane). 

In order to explain the long-term periodic variability of the ac- 
cretion rate, CG01 proposed a hierarchical triple stellar model, in 
which a distant tertiary exerts tidal forces on the inner binary. This 
results in a cyclic exchange of the angular momentum between the 
inner system and the tertiary (Kozai 1962; Lidov & Ziglin 1976; 
Mazeh & Shaham 1979; Bailyn & Grindlay 1987; Ford, Kozinsky 
& Rasio 2000, hereafter F00; Blaes, Lee & Socrates 2002; Wen 
2003, hereafter W03). The period of the resulting evolution of the 
parameters of the system is ~ P2/P1, where P2 is the orbital pe- 
riod of the tertiary, implying P2 ~ 1 d. The variable system param- 
eter relevant here is the eccentricity, ei, of the inner system, which 
causes changes of the distance between the inner Lagrange point, 
Li, and the center of mass of the secondary. This in turn changes 
the rate of the flow through L\ and the accretion rate. The mass of 
the tertiary has been constrained by CG01 as M3 < 0.5Mq based 
on the lack of its optical detection. We note that this constraint as- 
sumes the third body is a main-sequence star; if it is a white dwarf 
or a neutron star, M3 < 1.4Mq. (Still, the third star is most likely 
on the main sequence, which we assume hereafter.) However, no 
specific calculations of the proposed triple system, e.g., of the max- 
imum ei from modelling the variable rate of the flow through the 
Li point, or of the system parameters from modelling gravitational 
evolution of the system, have been done yet. 

Here, we first analyse the RXTE ASM monitoring data and the 
Proportional Counter Array (PCA) scanning data, and use them to 



quantify the source X-ray variability. We then present calculations 
on the dependence of the accretion rate through L\ on the eccen- 
tricity, which yields the maximum ei required to reproduce the ob- 
served amplitude of the superorbital variability. Then, we model 
evolution of hierarchical triple stellar systems, and reproduce the 
~170-d period and the maximum eccentricity implied by the L\- 
flow model. Our results put constraints on the masses of the system 
component, the inclination between the inner and outer orbits, and 
the ratio between its semimajor axes. We also report a discovery of 
a remarkable synchronization between the superorbital period and 
the period of the relativistic periastron precession of the binary. 



2 THE DATA 

Fig. [3 shows the long-term light curve of 4U 1820-303 from the 
RXTE ASM (1996 January 1-2007 February 20). For the clarity 
of display, we have averaged some adjacent 1-day measurements 
in order to achieve the minimum significance of 3a of the plotted 
count rate. We clearly see the cyclic variations of the count rate on 
the ~170-d quasi-period with a large relative amplitude. The ver- 
tical lines show the predicted minima according to the ephemeris 
of eq. (9) of CG01. Using the Lomb-Scargle periodogram (Lomb 
1976; Scargle 1982), we find that the present (daily-averaged) ASM 
data yield P3 = 170.6 ± 0.3 d, compatible with the periods of 
CG01 and ofW06. 

Fig.f2fa) shows the ASM light curve folded on the ephemeris 
of CG01, as well as the folded light curve averaged over 10 phase 
bins. We see that while the daily measurements span a factor > 10, 
the averaged light curve varies over the superorbital period span- 
ning a factor of ~2. 

We also use the available Galactic bulge scan dat|] from the 
RXTE PCA detector for the time interval of 1999 February 5-2006 
October 30. Fig. [2ft)) shows the count rate from those measure- 
ments folded, and folded over the superorbital ephemeris. We also 
show the count rate folded and averaged over 10 phase bins. Sim- 
ilarly to the ASM data, we see that while the individual measure- 
ments of the count rate span a factor ~10, the averaged light curve 
varies over the superorbital period spanning a factor ~2. We stress, 
however, that the variability is not strictly periodic, and P3 repre- 
sents only a quasi-period. Therefore, folding and averaging over 
a single period value results in some suppression of the actual 
average superorbital variability. With that caveat, we attribute the 
changes of the average flux within the factor ~2 to a relative sta- 
ble quasi-periodic process (Sections [3}0, and the remaining vari- 
ability to some aperiodic processes, in particular those operating in 
other atoll sources. We also assume that the ASM and PCA count 
rate variability reflects that of the accretion rate, and assume here- 
after its superorbital variability consists of variations within (0.7- 
1.4)(-M 2 ). 

Z07 have, in addition, analyzed available PCA and High En- 
ergy X-Ray Timing Experiment (HEXTE) pointed observations of 
4U 1820-303. Using a physically-motivated spectral model, they 
have estimated the bolometric flux from the source for each obser- 
vation. The resulting light curve, when folded over the ephemeris, 
is very similar to those of Fig. [2] In particular, the fractional vari- 
ability of the folded/averaged light curve is also found to be ~2. 
Since the bolometric flux is highly likely to be proportional to the 
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Figure 1. The RXTEI ASM light curve based on 1-day average measurements. The vertical lines show the minima of the superorbital cycle according to the 
ephemeris of CG01. Note departures of the ASM minima from the predicted dates up to ~50 d in some cases, as well as secondary minima (see also Simon 
2003). The solid curve presents our theoretical model of Sections|5]| 
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Figure 2. (a) The ASM light curve of Fig. [T] folded on the superorbital 
ephemeris of CG01. (b) The folded light curve for the PCA Galactic bulge 
scan data (see Section [2). On each panel, the histogram shows the light 
curve averaged over 10 phase bins, and the solid curve shows the theoretical 
model of Sections [5J2] (based on the second cycle of Fig. |4ji) fitted to the 
averaged data. 



rate of accretion onto a neutron star, that result confirms that the su- 
perorbital range of — M2 is indeed within a factor of ~2, as adopted 
above. 



3 THE MASS FLOW THROUGH THE INNER 
LAGRANGIAN POINT 

Here, we calculate the dependence of the rate of the Roche lobe 
overflow (assumed to equal the accretion rate) on the eccentricity. 
For that, we follow classical results on the Roche flow, expressing 
M2 as the product of the density and the sound speed at the in- 
ner Lagrangian point, L\, times an effective area of the flow. The 
position of the L\ varies with the orbital phase in an eccentric or- 
bit, yielding a varying M2. This yields the orbit-averaged M2 as a 
function of the eccentricity, ei. This in turn allows us to determine 
the maximum eccentricity, e max , required to enhance Mzie) by the 
factor of two inferred above from the varying Fboi, assuming that 
the minima of the superorbital cycle correspond to the e\ minimum 
of eo m 0. 

We first briefly estimate the parameters of the inner binary. 
We assume that the white dwarf fills its Roche lobe. Then, we com- 
bine the volume-averaged Roche-lobe radius for M2/M12 < 0.8 of 
Paczyhski (1971), 



-R2 _ _2_ ( M 2 \ 1/3 
ai ~ 3 4 / 3 \M 12 ) ' 

with the Kepler law, 

p2 

of = GM 12 —^, 
47I" 4 



(1) 



(2) 



where G is the gravitational constant, a\ is the semimajor axis, and 
M12 = Mi + M 2 . This yields the relation, 



Pi 



9tt R 



3/2 



2V2 (GM 2 ) 1 /2 



(3) 
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We then use the mass-radius relation for cold, low-mass stars of Za- 
polsky & Salpeter (1969; as fitted in R87), and assume pure He and 
the radius equal that of the fully degenerate configuration times a 
factor, / d . For / d = 1.1 ± 0.1 (R87), M 2 ~ O.O67lg;giEjM 
and7? 2 ^ 2.191q™ x 10 9 cm, where the lower and upper limits 
correspond to f c \ = 1 and 1.2, respectively. Hereafter, we assume 
/d = 1.1, for which the above M 2 and R 2 correspond to an He 
model of Deloye & Bildsten (2003) (who has calculated the mass- 
radius relation for white dwarfs of arbitrary degenerac}0) with the 
core temperature of ~ 10 7 K and the core density of ~ 10 4 ' 2 g 
cm -3 . We note that for a star approximated by an n — 3/2 poly- 
trope (Chandrasekhar 1939), we have 



A/2 



-1/3 



Rr. 



(4) 



r 2 a 0.0128(1 + xy /A — i 
V M 

(where X is the H mass fraction), which expression combined with 
equation 0} at X = yields M 2 ~ 0.067M Q and R 2 ~ 2.19 x 
10 9 cm, the values identical to those above for /d = 1.1. Then 
01 ~ 1.29 x 10 10 (Mi 2 /1.35Mq) 1/3 cm, where we used M 12 = 
1.35Mq corresponding to the presence of the resonance implied 
by equation l |38l > below, which then yields M\ ~ 1.28Mq. 

For those parameters, and assuming negligible effects of a 
possible outflow, eq. (17) of R87 yields the average mass trans- 
fer rate of (—Ma) ~ 3.5±o.g x 10 17 g s" 1 (where the lower and 
upper limits correspond to /d = 1, 1.2, respectively). This rate is 
fully compatible with the accretion rate, ~ (3.3±0.1) x 10 17 g s _1 
(where the error is statistical only), corresponding to the average 
bolometric flux measured using the pointed RXTE PCA/HEXTE 
observations, (F hol ) ~ (8.7 ± 0.2) x 10~ 9 erg cm" 2 s _1 (Z07) 
at D — 7.6 kpc and an accretion efficiency of 0.2. Using Mi = 
1.4Mq increases the above theoretical rate by ~ 0.2 x 10 17 g 
s" 1 . In our numerical estimates below, we use Mi = 1.28Mq, 
M 2 = 0.07Mq, R 2 = 2.2 x 10 9 cm, and the corresponding 
(-M2) = 4 x 10 17 gs' 1 . 

The distance, Rq, between the center of the white dwarf and 
the L\ point at e\ — 0, is given by the solution of 



/ -1 n-2 ,-, , M 2 _ 2 



(5) 



where x = Ro/ai, which yields R ~ 0.237ai ~ 3.05 x 10 9 
cm. In an elliptical orbit, the separation between the stars at a given 
orbital angle, <f> (measured with respect to the periastron), is given 

by, 

_ ai(l-e 2 ) 

s i - — r> (o) 



1 + ei cos <f) 
and the angular velocity is 

d(j> _ 2tt (1 + ei cos^) 2 
"dT ~~ A (1-e 2 ) 3 / 2 ' 



(7) 



At ei <C 1, the L\ distance is proportional to the separation, i.e., it 
equals Rq + d! , where 



d! ei + cos < 
K^- £l -- 



1 + ei cos ( 



(8) 



which changes between —ex and ei. (See Brown & Boyle 1984 for 
an expression valid at large ei.) 



2 We note that for a cold, fully degenerate, white dwarf with the parameters 
relevant to 4U 1820-303, the model of Deloye & Bildsten (2003) gives 
somewhat lower radii than that of Zapolsky & Salpeter (1969). 



The rate of the Roche lobe overflow can be generally written 
as (e.g., Savonije 1983), 



M 2 = Apc B , 



(9) 



where A is an effective area of the flow, p and c s are the mass 
density and the isothermal sound speed, c s = (kT/ prrm) 1 ^ 2 , re- 
spectively, at Li, and p is the mean molecular weight. Thus, M 2 
depends on the depth of the Li point within the donor star, d. It can 
be shown by expanding the effective gravitational potential around 
Li that 



A: 



Pt_GM 2 
2tt 



d. 



(10) 



(Savonije 1983). (Note that distinguishing here between the ra- 
dius of the L\ point, Ro, and the volume-averaged radius of the 
star filling the Roche lobe, R 2 , would require much more com- 
plex treatment of the flow than that adopted here.) Paczyhski & 
Sienkiewicz (1972) find that in the polytropic case —M 2 oc d 1,5+n , 
where n is the polytropic index [which can be shown to follow 
from equations j9TU0H. If we neglect for a while the illumination 
of the white dwarf by the X-ray source, the relevant index would be 
n = 3.25 of the (non-degenerate) surface layers of the white dwarf 
(Schwarzchild 1958). Then, the orbital-angle dependent accretion 
rate in an elliptical orbit can be written as, 



-M 2 (0,ei) 
M 



max(c(o — d' , 0) 



do 



(11) 



where do and Mo are the depth of the L\ point and the accretion 
rate, respectively, at ei =0 (and d = do — d'). Hereafter we 
assume that the minimum of the average light curve (see Fig. [2} 
corresponds to ei = 0, i.e., Mo = 0.7{—M 2 ) (see Section|2]l. 

For the assumed white dwarf parameters and the core tem- 
perature of 10 7 K (for the model of Deloye & Bildsten 2003; see 
above), we can calculate the intirinsic white dwarf luminosity [eq. 
(4. 1.1 1) of Shapiro & Teukolsky 1983]. For the metal abundance of 
Z ~ 0.01 of NGC 6624 given by Rich, Minniti & Liebert (1993), 
p — 4/3 for ionized He, and the (very approximate) Kramers' 
opacity of Schwarzschild (1958), we find that luminosity to be 
~ 10~ 3 Lq. Using the above values, we then use the photon diffu- 
sion equation, hydrostatic equlibrium, and the equation of state and 
solve for the profiles of T(d), oc d, and p(d), oc d 3,25 , of the non- 
degenerate white-dwar surface layer (Schwarzschild 1958; section 
4.1 of Shapiro & Teukolsky 1983). Then we find from equations 
l[9l00} that d ~ 1.64 x 10" 3 i? . At this d , T ~ 1.9 X 10 s K 
and p ~ 7.0 x 10~ 7 g cm -3 . Note that do is much less than the 
calculated thickness of the surface layer, ~ 0.01i?o, below which 
the electrons become partly degenerate at the interior temperature 
of ~ 10 7 K. 

The accretion rate averaged over the orbital period is, 



M 2 (ei) = 



(1-e 2 ) 3 / 2 
2tt 



M 2 



,ei) 



(12) 



(1 + ei cos (f>) 2 

Note that since this rate is averaged over time, a factor of dt/dcj), 
equation Q, appears in the integral over 4>. We will also use a di- 
mensionless accretion rate, 



-M 2 (ei) 
Mo 



(13) 



We find that a very low e max ~ 8 x 10 -4 is sufficient to increase 
—M 2 by a factor of 2. Fig. [3ja) shows —M 2 (<j), ei)/Afo of equa- 
tion i ll It at this value of ei. Fig.[3lb) shows the rh e of this model. 
Both the strong changes of the accretion rate over an orbit and the 
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Figure 3. (a) The dashed and solid curves show the orbital phase depen- 
dence of the accretion rate for the model without irradiation at ei = 
8 X 10 ~ 4 and do/Ro = 1.64 X 10 -3 , equation (TTJ, and with irradiation 
at ei = 4 X 10~ 3 and Ho/Ro = 2.3 X 10~ 3 , equation (16), respectively. 
The periastron corresponds to (f> = 0. In both cases, rh e ~ 2, illustrating 
the orbital phase dependence of Af 2 at the respective e max required to ex- 
plain the superorbital cycle of 4U 1820-303. See equation (7) f° r d<j)/dt 
(almost constant at ei <C 1). (b) The dashed and solid curves show the 
orbit-averaged accretion rate as function of e± for the two above models, 
respectively. 



fast increase of the orbit-averaged rate with ei are due to the very 
fast increase of the product pc B d with d, oc d 4 ' To , in the surface 
layer of the white dwarf (see above). We note that we have used 
here the structure of an isolated white dwarf, while the gravita- 
tional field close to L\ is affected by the gravity of the neutron star. 
This would flatten the density profile and lead to a requirement of 
a somewhat higher e max - 

Moreover, the white dwarf is very close to the X-ray source, 
and thus it is irradiated. Irradiation could be avoided if the accre- 
tion disk was strongly flared and obscured the white dwarf from the 
X-ray source. However, the observed strong orbital modulation in 
the UV (Anderson et al. 1997) is well explained by reprocessing of 
X-rays by the white dwarf (Arons & King 1993), which interpreta- 
tion rules out strong flaring. Thus, we consider now the case with 
irradiation of the white dwarf. 

The white dwarf subtends the solid angle ~ TvR^/af, which, 
for our binary parameters, is ~ 7.3 x 10 - 3 47r. For the (-Fbol) 
of Z07 (see above) and an albedo of 0.5 (e.g., Anderson 1981; 
London, McCray & Auer 1981), the white dwarf receives ~ 
2.2 x 10 35 (D/7.6 kpc) 2 erg s -1 , i.e., several orders of magni- 
tude more than its possible intrinsic luminosity. The absorbed flux 
per unit area at a mid-point of the white dwarf (at ~ ai — ro /2 



from the neutron star) varies then in the range Fi rl ~ (0.9- 
1.8) x 10 16 (D/7.6 kpc) 2 erg cm" 2 s" 1 . This variability is due 
to the X-ray luminosity undergoing the superorbital cycle with the 
amplitude of 2. The corresponding blackbody temperature changes 
from T = T ~ 1.1 x 10 5 (D/7.6 kpc) 1/2 K to 2 1/4 T . Assuming 
L oc m e , T(ei) = mi To. In optically-thick regions, the struc- 
ture of the atmosphere will be closely isothermal at this T (e.g., 
Anderson 1981). Thus, the atmosphere density will decrease ex- 
ponentially with the distance from the center of the white dwarf, 
p oc exp(— d! / H), where H is the scaleheight around the L\ point, 



H_ 
Ro~ 



<ZRo 
GM 2 



: 2.3 x 10~ 3 m e /4 



T 



1.1 x 10 5 K' 



(14) 



where we used our numerical values of the parameters for the sec- 
ond equality. Then, we use the results of Brown & Boyle (1984), 
who derived A — {2-k) 1 ^ 2 HRq (see their eq. 4), which can be 
expressed as, 



M 2 (0,ei) ~ (27r) 1/2 4 nC 3 p e X p(-d7H), 



(15) 



where tdyn = (R0/GM2) ' (— 55 s) is close to the dynamical 
time scale of the white dwarf [(7? 2 /GM 2 ) 1/2 ~ 30 s], and p is 
the atmosphere density at L\ corresponding to e\ — 0. Here we 
neglected the correction of Brown & Boyle (1984) for the velocity 
of the L\ point, since it is much slower than the sound speed in 
our case. Note that though the functional dependence oc cf po is 
most likely accurate, the constant of the proportionality remains 
relatively uncertain (see Savonije 1983) as well as the argument of 
the exponent may be more complex than that used above, see Frank, 
King & Raine (2002), pp. 352-353. The latter effect may increase 
the scaleheight around Li and thus lead to a requirement of a higher 
e ma x than than that estimated by us below. To solve this problem 
accurately, hydrodynamical simulations (see, e.g., Regos, Bailey & 
Mardling 2005) specific to 4U 1820-303 (beyond the scope of this 
work) are needed. 

The H and c s depend on ei through T(ei), and we denote 
their values at ei =0 as Ho and c s o, respectively. The accretion 
rate relative to that at ei = can then be written as, 



-M 2 (0,ei) 
Mo 



■ 3/8 

m e exp 



e\Ro ei 



H rhl /4 1 + ei cos < 



(16) 



This equation coupled with equations J121413t . giving m e , can be 
easily solved iteratively. 

This solution, for a given Ho / Ro, can be used to calculate the 
eccentricity required for a given amplitude of L. Fig.[3la) shows an 
example of the dependence of equation d 1 61 > for ei = 0.004 and 
Ho/Ro = 0.0023, yielding m e ~ 2 (i.e., e max ~ 0.004). Fig. 
Ob) shows -M 2 (ei) for the same Ho/Ro - These dependences are 
slower than those of the unirradiated case due to the slower (in 
the present case) increase of the density with the depth within the 
atmosphere as well its isothermality. 

In the irradiated case, the accretion rate enhancement due 
to the eccentricity depends primarily not on ei itself but on 
e\l (Ho/Ro). In particular, the eccentricity required to obtain the 
increase of the orbit-averaged accretion rate by 2 is e max ~ 
2Ho/Ro- Also, rh e increases initially, at low values of ei, very 
slowly. Thus, our assumption that Mo corresponds to ei = (see 
above) introduces only a slight error as long as e max 3> en. We 
also note that the value of e max does not depend on the (relatively 
uncertain) normalization of the — Af 2 dependence of equation l !15t . 

From that normalization, we can determine po as, 
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po 



Ah 



(2^/3^3, 



2.0 x 10~ 6 gcm~ 3 . 



(17) 



This corresponds to an electron density of no — 6.0 x 10 cm , 
and the Thomson optical depth above the L\ of 1.0. The Rosseland 
absorption opacity at the above p, T ~ 1.1 x 10 5 K and Z = 0.01 
is ~ 17 g -1 cm 2 (Iglesias & Rogers 1996), i.e., ~ 10 2 more than 
the scattering opacity for He. Thus, the medium is completely opti- 
cally thick, consistent with our assumption of the temperature equal 
to the blackbody one. The corresponding pressure is ~ 10 7 dyn 
cm -2 , which can be calculated to be ~ 10 2 times more than the 
critical pressure below which a transition to an optically thin regime 
begins (London et al. 1981). Note that the L\ region may be partly 
shadowed by the accretion disc, which would decrease the temper- 
ature and move the L\ region even more into the optically-thick 
regime. We also check that the time scale at which the atmosphere 
locally responds to the irradiation, ~ 2nokToHo / Fi tt ~ 15 ms, is 
much shorter than any other time scale of interest. 

On the other hand, the coefficient in equation J 1 5b and other 
details of the flow through L\ are relatively uncertain, and thus it 
is of interest to also consider the case with the L\ region being op- 
tically thin. Above a transition zone (Anderson 1981; London et 
al. 1981), the atmosphere temperature becomes close to the Comp- 
ton temperature of the surrounding radiation field (e.g., Kallman & 
McCray 1982; Begelman, McKee & Shields 1983), 

J huJ v Av 



(18) 



4k J J v dv' 

where J v is the mean intensity. For the sum of the X-ray spec- 
trum of the source, its reflection from the star, and the blackbody 
emission from reprocessing in the underlying regions of the star, 
Tq will be ~ 2/3 of the Compton temperature for the X-ray spec- 
trum alone. Based on the spectral fits of Z07, we have calculated 
the latter to range from (1.6-1.8) x 10 7 K in the high-luminosity 
state to ~ 10 8 K in the low-luminosity state. Thus, Tc ~ 10 7 
K for the dominant high-luminosity state. As the Compton tem- 
perature depends only on the shape of the spectrum but not on its 
flux, T ~ Tc is almost independent of — Ah, except for the low- 
est — M2 corresponding to the hard low state, where Tc is several 
times higher. Neglecting the last effect (important only at the flux 
minima), equation d 14b yields H / Rq ~ 0.2. 

If this condition were dominant over the superorbital cycle, a 
rather high e max ~ 0.4 would be required to account for the su- 
perorbital variability. This would require changing the radius (and 
mass) of the white dwarf, and raise the issue of the stability of the 
mass transfer (beyond the scope of this work). However, we stress 
that according to our estimates the L\ region of the irradiated star 
is unlikely to be in the optically thin regime. 

On the other hand, we cannot rule out that the L\ point is in 
the transition zone (Anderson 1981; London et al. 1981) between 
T ~ 10 5 K and 10 7 K, in which case detailed solutions of the atmo- 
sphere structure would be required to calculate — Mz{e\), and the 
e m ax would be somewhat higher than the value of 0.004 estimated 
above. On the other hand, the shadowing of the white dwarf by the 
disc would (as mentioned above) decrease the temperature of the 
L\ region. These effects introduce some systematic uncertainties 
on the value of e max . 



4 MODELS OF THE TRIPLE SYSTEM 

Above, we have modelled the accretion rate variability in 4U 1820- 
303 as due to variability of the binary eccentricity. Here, we model 



the required variability of the eccentricity as due to gravity of a 
third star in the system. The third star should have the semimajor 
axis (0,2 ; measured with respect to the center of mass of the inner 
binary) satisfying 02/01 S> 1 in order not to perturb the binary 
motion on time scales in the range between Pi and P3, which per- 
turbations have not been observed. Thus, this model is of a hier- 
achical triple system. We then attempt to reproduce the observed 
171-d period and its relative coherence (Section |2]l as well as the 
eccentricity amplitude necessary to reproduce the observed ampli- 
tude of the superorbital variability as inferred from our Li-flow 
calculations (Section[3j. 

We begin with a brief review of the main features of secular 
effects in hierarchical triple systems. It has been known (e.g., Kozai 
1962) that a third outer body can induce quasi-periodic oscillations 
of the inner eccentricity with a long quasiperiod of 



K 



Pi 



where P2 (the outer orbital period) is given by 



,2 47r 2 a| 



(19) 



(20) 



" 2 ~~GM~' 

M = Mi + Ai 2 + M 3 , and K (often ~1) depends on the sys- 
tem parameters. There are two main regimes. At an initial mutual 
inclination, io, above a critical angle, i c — 40°, the amplitude of 
the inner eccentricity is large and K ~ 1. The evolution can be 
well described by taking into account only the lowest-order per- 
turbative term in the system Hamiltonian expanded in powers of 
ffll/02, namely the quadrupole term, oc (ai/a2) 2 . Then, the mu- 
tual inclination, i, is anticorrelated with the inner eccentricity, ei, 
with (1 — e 2 ) cos 2 i ~ constant. In this approximation, both semi- 
major axes, the outer eccentricity, e%, and the magnitude of the 
angular momentum of the outer binary are constant, and there is 
exchange of the angular momentum between the inner and outer 
binary. In a conservative system, also the total angular momen- 
tum and energy are constant. The maximum inner eccentricity is 
rougly e max ~ [1 — (5/3) cos 2 io] 1 ^ 2 (Innanen et al. 1997; Hol- 
man, Touma & Tremaine 1997) above the critical angle (and thus 
e max ~lfori ~90°). 

Below i c , e max <C 1, and the quadrupole approximation be- 
comes insufficient. In particular, for coplanar orbits, io = 0, the 
quadrupole term in the Hamiltonian becomes null, and the first non- 
zero term is octupole, oc (ai/a2) 3 . In the octupole approximation, 
only the semimajor axes are constant (apart from the total energy 
and angular momentum). In this regime, the inner eccentricity also 
varies quasi-periodically, but at a period longer than that of the 
quadrupole one. In the intermediate regime where both terms are 
important, time dependencies of the inner eccentricity show both 
periodicities (Krymolowski & Mazeh 1999). Also, at e max <C 1 in 
general, e max decreases with the increasing 0,2/0,1, increases with 
the increasing the initial outer eccentricity, e2,o, and it depends very 
weakly on either M 3 /Mi 2 or M 2 /Mi (F00). 

We use a numerical model calculating the time evolution of an 
isolated hierarchical triple of point masses, using secular perturba- 
tion theory up to the octupole terms (e.g., F00; Miller & Hamilton 
2002; Blaes et al. 2002; W03). We neglect possible effect of the 
tidal deformation of the white dwarf on the gravitational force ex- 
erted on the inner binary, see, e.g., Soderhjelm (1984), Eggleton & 
Kiseleva-Eggleton (2001). Also, the Hamiltonian is averaged over 
the periods of both the inner and outer binaries, and thus any pos- 
sible short-time scale changes (e.g., Bailyn 1987; Georgakarakos 
2002) are averaged out. Our calculations are Newtonian apart from 
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the general relativistic (GR) periastron precession in the inner bi- 
nary. Its first-order post-Newtonian period, Ppn, is given by (e.g., 
Weinberg 1973), 



Ppn = 



P iai c 2 (l-el) _ P^VCl-e?) 



(21) 



3GMi2 3(27tGMi 2 ) 2 / 3 ' 

We include this effect in the same way as in, e.g., W03. 

The free parameters of the model are the masses of the neutron 
star, Mi, and of the third star, M3, 02/01, io, ea.Os and eo- We 
assume 0<eo >C e max ~ 0.004 (Section[3j. Then, we know Pi 
and P 3 , and assume M 2 = 0.07Mq (Section[3}. 

Given the relatively large number of free parameters and com- 
plex dependencies between them and the resulting period and the 
amplitude of the eccentricity modulation, it has proven difficult 
to constrain the parameter space numerically. Thus, we have ob- 
tained some approximate analytical constraints. Given the presence 
of only one, well-defined, long-term periodicity in the system (Sec- 
tion [2j W06), evolution of the triple system can be dominated by 
either the quadrupole or the octupole term, but not by both, which 
would have given rise to two periodicities (as noted above). Given 
the considerably greater simplicity of the quadrupole equations, we 
consider them for our analytical estimates (though we do include 
the octupole term in numerical calculations). Given our require- 
ment of e max <C 1, this implies values of io close to i c . 

We first estimate the maximum eccentricity. We use the con- 
servation of the total angular momentum and energy. The former 
gives us the mutual inclination, i, in terms of io and ei, e.g., eq. (3) 
in Miller & Hamilton (2002). Then the (quadrupole) Hamiltonian 
gives the total energy. The minimum and maximum inner eccen- 
tricity corresponds to the inner periastron angle (measured from 
the intersection of the two orbits), gi, of and 7r/2, respectively. 
Equating the Hamiltonian at those two angles and using i from the 
angular momentum conservation, gives, in the lowest (second) or- 
der of ei, 



4 + 6W + 2cosi //3 
10 cos 2 io — 6 + #pn + 2 cos io / f3 



2 



where 



8GM? 2 aZ(l - e 2 2 



\3/2 



MU 2 M 3 a^(l-ei ) 



c 2 M 3 ai 

1/2/ 



\l/2 



M 1 / 2 M 1 M 2 a 1 



1/2 



(22) 



(23) 



(24) 



In order for a solution with e max ^ eo to exist, the denominator 
in equation J22b should be nearly zero. Also, /3 > 1 at 02 > di, 
which we assume hereafter. Thus, 



2 • 2 • 

COS to — COS l c 



6 - 6> F 
10 



(25) 



Note that this i c is also the critical angle for the large 6max regime 
(see Blaes et al. 2002), and it becomes the Newtonian critical angle 
when #pn = (Kozai 1962). In the small e\ limit, the initial mu- 
tual inclination angle should be approximately equal to but slightly 
smaller than the critical value, which is different from the require- 
ment for high e max regime. The value of e max is determined by 
how close i is to i c . 

For a solution of equation ( 125b to exist, (9pn < 6 is required. 
Then, we have from equation ( 1231 ), 



al < 3 c 2 aiM 3 



a? ~ 4 4GM?- 



1 



„2 N-3/2 
e 2,0) 



(26) 



For M12 ~ 1.5, the upper limit on M 3 of O.5M (CG01) and 



e| <S 1, 02/01 < 25. Note that the constraint of equation l !26t is the 
same as that derived for the high-ei case [except for the (1 — e 2 ) 3 ^ 2 
factor, Blaes et al. 2002]. It is related to the fact that the Newtonian 
regime corresponds to Ppn 3> P3. Otherwise, the GR periastron 
precession decreases the amplitude of the Kozai oscillations of ei, 
and leads to its disappearing in opposite limit. This is because the 
effect of the third body is a coherent summation of the tidal per- 
turbations over many orbital periods, and the GR precession partly 
destroys this coherence (e.g., W03). In particular, this leads to a 
change of i c , equation J25t . reducing the high-e max regime. 

Another constraint on the parameter space comes from the 
observed P3 ~ 171 d. The inner binary spends most of the 
time around gi = 7r/2, at which dgi/dt — 0. We thus write 
gi — tv/2 + 8 with 5 <C 1 and expand the (quadrupole) evolu- 
tion equations for ei and g\ [e.g., eqs. (16-17) in W03] to the first 
order in 8 and in the small-ei approximation. We then use equa- 
tion i22l to express the results in terms of e max /eo, divide dei/dt 
by dgi/dt, and integrate over e\ from eo to e max . This yields the 
value of K [equation l |19H of 

2 5/2 / M 

, , / . . ^ — r(i -c 2 n; - 

eo 



K 



2 \3/2e m ax , 1 e max 
2.nJ , 



(27) 



3tt M 3 (4 + pn ) ^" e e 

where / is a fudge factor to compensate for the approximation we 
used in the derivation. We found empirically that / — 1 are within 
a factor of two, and in our example described below, / ~ 1.1. Note 
also that the superorbital period is strongly dependent on eo, on 
which the accretion rate depends very weakly (as long as eo — 0, 
Fig.[3p). 

The following constraint on the total mass of the inner binary 
can be obtained using equations l |19l420t , l |23t and ( |27t , 

3(GM /c 3 ) 2 / 3 (2^) 5 / 3 P 3 e o _ /ff PN / Afa \ ~ 2/z 
2V2 Pl 5 / 3 emax l n V2 (emax/eo) 4 + PN ^M ; ' 

which (for the observed Pi and P3) can be solved for #pn as, 

4 



#PN 



0.274/ (M12/M0) 



-2/3 



• In 2 



1 



(29) 



The constraint of #pn < 6 then yields a relation between M12 and 

,: max /eo, 

3/2 

(30) 



§1 < 0.067 



ln 3/4 



eo 



At / ~ 1.1, it allows M12 1.5M provided e ma x/eo ^ 5.5. 
Since the required relative amplitude of AI2 of ~ 2 can be achieved 
at any e max /eo > 3 (see the solid curve in Fig. [3J5), this is only a 
very weak constraint, allowing practically any of the theoretically 
possible masses of the neutron star at modest values of e max /eo. 

Equation d23t at #pn < 6 also yields a constraint on M3. 
If we assume M12 — 1.5Mq and that the system is hierarchical, 
a 2 /ai> 5, we obtain M3 > 0.004Mq [see also equation 126H , An 
independent relation follows from equations dl91420t and ( I27t . 



M 3 



2 5 / 2 a 3 Pie n 



(l-e 2 )3/ 2 Mi2 37ra?(4 + 0p N )P 3 eo 



ln*^ 



eo 



(3D 



which yields the same constraint of M3 > 0.004Mq at 
e ma x/eo > 7. Thus, even a very low-mass third body can 
still induce the required eccentricity oscillations in the inner 
binary. 

We have not studied analytically constraints in the octupole- 
dominated regime (which regime has been considered, e.g., by Ra- 
sio 1994, 1995; F00; Georgakarakos 2002; Lee & Peale 2003). The 
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octupole-dominated regime may likely yield solutions with io ~ 
(as compared to io — i c in the quadrupole regime). Still, if such 
solutions exist, our conclusion of only very weak constraints on 
the masses of the system components would remain unaffected. 
On the other hand, an important difference between solutions in 
the two regimes is the variability of the mutual inclination. In the 
quadrupole low-ei regime, the maximum change, Ai — i — io, is 
given by 



Ai 



2 2 



cot io, 



(32) 



which implies only a very small change of i over the course of the 
superorbital cycle, e.g., about —2.3" at io = 40°, e max = 0.004. 
On the other hand, numerical results of F00 (fig. 8) show Ai ~ 
±15° in their example for the octupole-regime, low-ei oscillations, 
for which io = 15°, e max ~ 0.001. Since the angular momentum 
is dominated by the outer binary, i is approximately equal to the 
inclination with respect to the constant axis of the total angular 
momentum, which, in turn, is related to the value of (in principle 
observable) the inclination of the inner orbit to the line of sight. 

For assumed values of Mi, M2, M3, eo, e max , and with the 
observed Pi and P3, the above equations can be used to find a-i 
and io- These appproximate analytical solutions can be then fine- 
tuned numerically. We first present an example yielding P3 ~ 171 
d within ~1 per cent and e max — 0.004, corresponding to our 
model of the accretion with irradiation (Section [3}. Its parameters 
are Mi = 1.29M , M 2 = O.O7M , e = 10" 4 , M 3 = O.5M , 
a 2 /ai = 8.66 (yielding P 2 ~ 0.17 d), i = 40.96°, e 2 ,o = 10" 4 , 
corresponding to #pn — 0.22, j3 ~ 19, and K ~ 41. The initial 
values of the periaxis angles, gi and g 2 , are set to zero. Note that 
e max /eo is very sensitive to io, see equation l !22t . 

Time evolution of this model is shown in Fig. |4ja). Fig. [5] 
shows the Fourier power spectrum of ei at the sampling rate of 
20 d _1 (to approximate that of the RXTEIhSM). The P 3 peak is 
unresolved at the resolution of the Fourier spectrum. This narrow- 
ness is similar to that observed, cf. fig. 19 of W06. In Fig.[4fa), we 
also see a weak second quasi-periodicity, with the period slightly 
longer than twice of the main one. This is an effect of the octupole 
term in the evolution equations (see above). 

We then use our results of Section [3] connecting ei to — M2. 
We first apply equations dl6| > and d!2H13t . yielding — Ma(e) with 
the flow parameters as in Fig. [3] to convert ei(£) of Fig. |4fa) to 
—Mi(t). We multiply the model period by 0.987 in order to obtain 
the exact agreement with the observed P3 . We first plot the result- 
ing — Miif) (fitting only the normalization) by the solid curve in 
Fig.Q] We see that the model fits well the overall superorbital vari- 
ability. We also note that the exact Newtonian calculations of the 
evolution would increase the scatter in ei(t), i.e., leading to a less 
periodic behaviour (e.g., F00). This could explain the presence of 
some observed superorbital eyeless with a significantly different 
duration that the average P3 (see Fig.QJ. 

We then take the second cycle in Fig. |4{a) as representative 
and fit it to the average superorbital phase diagrams of Figs.|2ja-b), 
with the second free parameter being the phase offset (with respect 
to the ephemeris of CG01), which we find <^3,o/2-7r = —0.088 and 
-0.076 for the ASM and PCA data, respectively. The fit results, 
shown by the solid curves in Fig.[2la-b), reproduce well the aver- 
age phase diagrams. 

Finally, we also find a triple solution corresponding to our 
model without irradiation of the white dwarf, which requires 
e max ~ 8 x 10~ 4 (Section [3}- Its parameters are identical to the 
previous model except for 02/01 = 16.95 (yielding P2 ~ 0.47 
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Figure 4. The time evolution of the inner eccentricity, e\ , for our best model 
(a) with irradiation of the white dwarf and (b) without irradiation. 




Figure 5. The Fourier power spectrum of the eccentricity dependence for 
the case with irradiation, Fig.|4}a). 



d) and io = 48.5°, corresponding to #pn — 1.65, j3 — 26.5, and 
K ~ 4.6. Time evolution of this model is shown in Fig.|4jb). 



5 LONG-TERM EVOLUTION AND OTHER 
CONSTRAINTS 

Above, we have neglected the effects of emission of gravitational 
radiation (and of any other process of long-term loss of angular 
momentum) on the evolution of the inner binary. This is fully jus- 
tified as we have considered evolution over the course of a decade, 
which is many orders of magnitude shorter than the time scale due 
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to emission of gravitational waves. Thus, we can consider the long 
term evolution separately from the short-term one. 

The time scale for the loss of the angular momentum, Ji , from 
the inner system due to emission of gravitational radiation is given 
by (Peters 1964), 



is. 

Ji 



32 
5 



27r y/3 G 5/3 MiMa 1 + ( 7 /g) e 2 

c 5 MU 3 (l~e?) 5/2 



1.23 x 10" 



Pi. 
MiM 2 



(33) 



M^ 3 M 



1/3 
12 



-8/3 



where 



Ji = MiM 2 



Gai(l-e?) 



Mi; 



1/2 



ei < 1, (34) 



(35) 



This yields r g ~ 1.1 x 10 yr for our system parameters. We see 
that the effect of ei > on the rate of emission of gravitational 
radiation is negligible as long as e\ <C 1, 

By logarithmically differentiating equation {35}, assuming 
conservative mass transfer (Mi = — M 2 ) and using equations {TJ- 
[2} and {4}, we can relate Ji to M 2 and A (assuming e\ <C 1), 



(36) 



Jl 


-M 2 




M 2 \ 


Pi 




M 2 \ 


Jl 


M 2 


(I- 


Mi J 


= PT 


(§- 


Mi J 



where J\ includes all contributions to the inner angular momentum 
loss. Note that although the white dwarf expands as it loses mass, 
the accretion from M% is not self-sustaining, i.e., M 2 = if J = 0, 
as well as it can take place only for M 2 /Mi < 2/3. Also, Pi /Pi = 
— M 2 /M 2 > 0. This is in conflict with the observations, showing 
Pi < (see Section [T}. We note here that taking into account a 
possible mass loss from the system cannot resolve this discrepancy. 
If a fraction, 5, of the mass lost by the secondary is ejected from the 
system, we derive from equations {TJ{2} and {4}, 



P 
Pi 



-M 2 
M 2 



■ M 2 
'Mii 



(37) 



Thus, although the mass loss reduces Pi, it is still > 0. Only in 
the absence of mass transfer, when the companion does not fill its 
Roche lobe, angular momentum losses (in particular, those due to 
gravitational radiation) lead to Pi /Pi — 3 J / J < 0. A possible 
resolution of this conflict is the apparent observed Pi < resulting 
from gravitational acceleration of the system in the cluster potential 
(CG01). 

Then, the likely evolutionary scenario of 4U 1 820-303 is still 
that of R87, as shown on their fig. 1. According to it, the period 
shortly after the onset of the mass transfer (~ 10 6 yr ago) was ~500 
s. Then, oi was lower and #pn higher [equation [23}] than now. 
Consequently, cos 2 i c was less than the present value (i.e., the low- 
e ma x, octupole, regime was increased), see equation ((25}. Thus, the 
eccentricity oscillations had in the past lower amplitudes than now, 
and the system was in the octupole-dominated regime. On the other 
hand, a future increase of Pi may lead to i > i c and moving the 
system to the high-e max regime (provided io ~ i c at present). 

On the other hand, we find that at the present moment of the 
source evolution, the GR periastron precession period, Ppn, is very 
close, and may be equal exactly, to the superorbital period, P3. 
Namely, 

5 / 3 / Mv N " 2/3 



PPN „ .2, / Pi 



(i-O 



-V 



(38) 



171 d v ""'\685sJ V L347M 0/ 
We consider this equality to be a very remarkable coincidence. We 



do not see any straightforward way in which the GR periastron pre- 
cession could by itself (i.e., without the presence of a third star) af- 
fect the orbit-averaged accretion rate. Thus, the origin of the above 
equality can be either purely accidental or it may be due to some 
evolutionary processes not understood yet. We note that F00 find 
the presence of a resonance at Ppn — Pi manifesting itself as a 
peak in e max (fig. 14 in F00). However, it is not clear how that 
resonance could lead to the synchronization of Ppn — P3. 

On the other hand, our results and those of R87 do allow 
Ppn = P3 to be satisfied exactly. If this is the case and at 
M 2 ~ (O.O6-O.O8)M , ei < 1 (Section[3}, the mass of the neu- 
tron star would be Mi ~ (1. 27-1. 29)Mq. Remarkably, modelling 
of an X-ray burst of 4U 1820-303 by Shaposhnikov & Titarchuk 
(2004) gives Mi = 1.29lg;^M . Also, Grindlay et al. (1984) 
have measured the total masses of low-mass X-ray binaries in glob- 
ular clusters (including 4U 1 820-303) and concluded that the ini- 
tial neutron-star masses appear substantially less than 1.4Mq on 
average. We also note that some low neutron-star masses have re- 
cently been measured with high precision, e.g., 1.250+q '005 Mq 
for PSR J0737-3039B (Lyne et al. 2004), and 1.18±^M for 
the neutron-star companion of PSR J1756-2251 (Faulkner et al. 
2005). Given those results, Mi = 1.28Mq in 4U 1820-303 ap- 
pears entirely plausible. 

We also note that although Arons & King (1993) proposed that 
Mi ~ 2Mq if the initial mass of the white dwarf were ~ 0.6Mq 
and ~ 0.5Mq has been accreted, the calculations of R87 imply 
a much lower initial mass and a short mass-transfer time interval, 
with only <C O.IMq accreted. Then, Zhang et al. (1998) proposed 
that Mi ~ 2.2Mq if the frequency of the upper kHz QPO at its ap- 
parent saturation at 1060 Hz were equal to the Keplerian frequency 
at the last stable orbit. However, van der Klis (2000) shows that 
there is no saturation in a more extensive data set. Also, current 
theoretical models of kHz QPOs do not postulate that frequency 
identification, and use it only to provide an upper mass limit (e.g., 
van der Klis 2000). 

The presence of the third star may affect the system evo- 
lution only if it remains stable and is not disrupted by encoun- 
ters with stars in the cluster. The triple system itself is stable if 
a 2 /oi > 2.8(1 + M 3 /Mi 2 ) 2/5 (Mardling & Aarseth 2001), which 
is clearly satisfied in our case, with a 2 /ffli ~ 10, M12 ~ 1.5Mq 
(Section!), and M 3 <0.5Mq (CG01). Then, the encounter time 
scale has been estimated by Miller & Hamilton (2002) in the limit 
of domination of gravitational focusing, and we write it as 



: 4 x 10 



8 4 x 10 5 pc" 3 lO n cm2M 

IlGC 02 M 



■yr- 



(39) 



Here, ugc is the number density of stars, which we assume to be 
numerically equal to the mass density of ~ 4 x 10 5 Mq pc~ 3 
in the core of NGC 6624 (Ivanova et al. 2005), and a 2 ~ 10 11 
cm, M ~ 2Mq (Section |4}- Thus, r onc is much longer than the 
gravitational-radiation time scale, equation {34}, but still shorter 
than an estimated cooling time of the white dwarf in 4U 1820-303 
of ~ 10 9 yr (the appendix of R87). Thus, an encounter could have 
formed the triple system during the lifetime of 4U 1820-303, but it 
is stable on the current evolutionary time scale. 

We also note that the angular momentum will also be lost from 
the system due to tidal dissipation within the secondary at ei > 0. 
This mechanism usually leads to circularization (and synchroniza- 
tion) of close binaries (e.g., Zahn 1977), but in our case the in- 
fluence of the tertiary will force (ei) > and continuous dissi- 
pation (Mazeh & Shaham 1979). The time scale for this process, 
Tt = —Ji/Jt, is related to the circularization time scale. This, 
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unfortunately, remains very uncertain for stars without convective 
envelopes, in particular white dwarfs (e.g., Zahn 1977), and it may 
be very long. Here, we consider the turbulent dissipation model of 
Press, Wiita & Smarr (1975), which probably gives the the lower 
limit to the actual dissipation time scale in a white dwarf (Zahn 
1977). Using the time scale of Press et al. (1975) and the approach 
of Mazeh & Shaham (1979), we obtain, 



n 



125tt (1 - e( 
121 ~~ 



2-,ll/2 



a.1 \ Rt 



Rz J K„ M*M 12 



Pi, 



(40) 



where Rt is the effective Reynolds number, which we assume = 
20 following Press et al. (1975), and is a dimensionless factor 
of the mean turbulent viscosity, 



Ku = 



14 



M2RI 1 



p(r)r 13 dr, 



(41) 



where p(r) is the stellar density. We have calculated ~ 0.0257 
for a polytrope with n = 3/2, appropriate for a low-mass white 
dwarf. For our system parameters, we find then r t ~ 2300/ (e 3 ) 
yr (where the average is over the superorbital cycle), which for 
(e?) 1/3 = 0.003 becomes ~ 10 11 yr. Thus, even if the dissipation 
is as fast as envisaged by Press et al. (1975), the process is com- 
pletely negligible compared to the gravitational radiation. However, 
n <r g for (e?) 1/3 > 0.06. 

On the other hand, the soft X-ray absorption towards the 
source is consistent with being caused entirely by the interstellar 
medium, with no evidence for any local component (Futamoto et 
al. 2004). This appears to rule out strong outflows from either the 
irradiated disc, the irradiated surface of the white dwarf, or, in par- 
ticular, mass loss through the L2 point. At the mass ratio of the 
binary, the absence of the L2 mass loss by the Roche-lobe filling 
star imposes the constraint of e max < 0.07, as calculated by Regos, 
Bailey & Mardling (2005). This rules out the above condition of 
(e 3 ) 1 ^ 3 > 0.06. Thus, presently, tidal dissipation appears not to be 
the dominant channel of angular momentum loss. 



6 CONCLUSIONS 

We have obtained the following main results. 

1 . Using the RXTE data, we quantify the average amplitude of 
the superorbital variability of the source. We find that amplitude in 
the folded and averaged light curves to be by about a factor of two. 

2. We consider the model in which the superorbital variability 
is due to a third star inducing cyclic variations (Kozai 1962) of the 
eccentricity of the binary (as proposed by CG01). We have stud- 
ied the dependence of the rate of the Roche-lobe overflow on the 
eccentricity taking into account the strong irradiation of the white 
dwarf by the X-ray source (Fig.[3p). We find the amplitude of the 
eccentricity required to account for the variability of the accretion 
rate by the factor of two is e max ~ 0.004 (assuming the minimum 
of the eccentricity is close to null). However, systematic uncertain- 
ties of this model may somewhat change the actual value of e max . 

3. We reproduce the above maximum eccentricity of ~ 0.004 
and the observed superorbital period of 171 d in a detailed model of 
a hierarchical triple system. Our calculations are Newtonian except 
for inclusion of the GR periastron precession. We then convolve the 
obtained eccentricity light curve with our theoretical dependence 
of the accretion rate on the eccentricity and obtain a theoretical 
luminosity light curve. This theoretical light curve is compared to 
and found to be in a good agreement with the observed light curves 
from the ASM and PCA detectors, see Figs.[[H2] 



4. We also study the parameter space allowed by the obser- 
vational data, the e max ~ 0.004 constraint, and A/3 < 0.5Mq 
(CG01). We obtain analytical solutions for the low-eccentricity 
Kozai oscillations in the regime dominated by the quadrupole terms 
of the evolution equations. We find the ratio of the semimajor axes 
to be < 25, which follows from the constraint that the GR peri- 
astron precession does not quench the Kozai eccentricity modula- 
tion. The masses of the system components are only weakly con- 
strained. In particular, the canonical neutron-star mass of 1.4Mq 
is allowed in this model. The mass of the third body is constrainted 
as 0.004 <Af 3 /M Q < 0.5. 

5. We find that the period of the binary GR periastron preces- 
sion is approximately equal (and it is allowed to be exactly equal) 
to the observed superorbital period (which equals the period of the 
Kozai eccentricity oscillations in our model). We find this to be a 
remarkable and puzzling example of synchronization in a physical 
system. 

6. We obtain some other constraints on the system. The bi- 
nary eccentricity has to be < 0.07 from the apparent absence of 
the mass loss by the Roche-lobe filling white dwarf through the 
L2 point. The angular momentum loss due to tidal dissipation in 
the white dwarf is found to be negligible compared to the loss due 
to emission of gravitational radiation. Also, we find the theoretical 
mass transfer rate due to the angular momentum loss via gravita- 
tional radiation is in complete agreement with that corresponding 
to the average bolometric flux from this source (as measured by 
Z07). 
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